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Abstract 

We present a computational random walk method of measuring the isother- 
mal pressure of the lattice gas with and without the excluded volume interac- 
tion. The method is based on the discretization of the exact thermodynamic 
relation for the pressure. The simulation results are in excellent agreement 
with the theoretical predictions. 

I. Introduction 

The pressure is defined as the force per unit area and is one of the measurable macroscopic 
thermodynamic quantities. It is a macroscopic manifestation of the microscopic phenomena 
of atomic collisions, by which atoms or molecules transfer the momentum and thus exert 
the force to the wall. At the macroscopic level, the pressure can be directly measured by 
monitoring or recording the force on the wall of a container. At the microscopic level, 
a conventional way is to design a method that keeps track of the momentum transfered 
by atoms to the wall via collisions, and this may be done by Molecular Dynamics(MD) 
simulations.^ However, developing and learning MD simulations is time consuming, and 
is often well beyond the scope of undergraduate or introductory graduate level statistical 
mechanics. To the best of the authors' knowledge, no conventional introductory statistical 
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mechanics text books have yet dealt with the computational methods of measuring the 
pressure. Considering the surge of computer-assisted course developments in recent years, it 
may be desirable to design a computational method of measuring the pressure. The purpose 
of this paper is to introduce such a Monte Carlo method, which is conceptually simple and 
extremely easy to implement with a minimal knowledge of comutation, because it only deals 
with the random walk of point particles on a lattice. We first explain the mathematical 
basis of the method in section II and III, and then present simulation results for the lattice 
gas with and without excluded volume interaction in section IV. 

II. Background 

Our statring point is the mathematical expression of the pressure P of a container of 
volume V at a given temperature T, which can be found in any standard text book of 
statistical mechanics^, 

P = -dF/dV = (1) 

oV 

where F — —kTlnZ is the free energy of the system with Z the partition function and k the 
Boltzmann constant. The contribution to the partition function comes from two sources, one 
from the kinetic energy and the other from the potential energy. The latter is known as the 
configurational integral, which invovles the integration of interaction energy among particles. 
The former is a function of temperature, and is independent of the volume. Therefore, 
at equilibrium, the contribution coming from the kinetics is decoupled, and is separated 
out from the configurational statistics, and the equilibrium pressure can be determined 
by considering only the configurational properties of the system. The nature of interactions 
among particles is determiend by the form of the potential energy. In this paper, we consider 
the interaction potential U to be either nonexistent, in which case the pressure is given by the 
ideal gas law, or the particles interact with each other via hard sphere potential. Even with 
hard sphere potential, the exact expression of the pressure for the continuum system is not 
yet known, even though a few approximate closed forms for the pressure have been reported 
in the literature. ^'^'^ However, for the lattice gas, the situation becomes much simpler. 
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because in this case, the free energy can be computed exactly. In the lattice version, the 
hard sphere interaction may be equivalent to introducing the excluded volume interaction 
by prohibiting the multiple occupancy of particles on a given lattice site. Without such 
restriction, the system will follow the ideal gas statistics. Note that in both cases, the 
contribution to the free energy comes only from entropic part because the interaction energy 
is zero. Hence, F = —TS = —kTlnZ, with S the entropy of the system. We now compute 
the entropy of the lattice gas, the result of which will then be used to test the validity of 
the random walk method. 

Consider now a lattice model, where N particles are confined to move only on d dimen- 
sional discrete lattice points of a container with the volume V — H • L^~^ with H the height 
and L*^""^ the area of the side wall, where the pressure is being monitored. Note that we set 
the lattice spacing to the unit value. In the absence of excluded volume interactions, the 
entropy of the system is given by: 

S = klnVL = kln{uj^) (2) 

where klncu is the entropy of a single particle, which is simply the volume V. Hence, we 
obtain the ideal gas law. 

P ^ -T{^) ^ kyT ^ k<f>T (3) 

with 4> = N/V the density of the particle. 

In the presence of the excluded volume interaction, the entropy of the system is no longer 
the product of one particle entropy. Rather, it is the total number of ways of putting N 
particles in a volume V. In the discrete version, it becomes, 

S = klnVLiV, N) ^ -V[(t)ln(t) + (1 - 0)/n(l - 0)] (4) 

where Q(y, N) = V\/N\{V — N)\ and the Sterhng's formula has been invoked. The equilib- 
rium pressure is then given by: 
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kTln{l - (j)) 



(5) 



One may use the grand partition function approach to obtain the similar result^. In this case, 
the grand partition function, Q, becomes: Q = I]^=o = (1 + 2;)^, where the fugacity 

z = exp{j3^) and the canonical partition function = V\/N\{y — N)\. Now, from the 
relations (see Chapter 7 of Huang in ref.2) PV/kT — InQ, and N — zdlogQ/dz — zVj (1+-^), 



the logarithmic singularity at the closed packed density. For the continuum system, it has 
been proved that the pressure exhibits the algebraic singularity.^ 

III. The Mathematical Basis for a Random Walk Monte Carlo Technique 

We now describe a mathematical basis for a Monte Carlo technique that computes the 
pressure of the lattice gas numerically. If the method is correct, then it should reproduce Eqs. 
(3) and (5). This is based on a method derived first for the continuum system by Percus^, and 
then extended to a lattice system to measure the pressure in polymeric fluids simulations.^ 
For the sake of completeness, we will briefly describe the essence of the formalism here, 
but the reader is referred to ref.7 and especially ref. 8 where the essential idea has been 
advanced. For a d dimensional lattice gas in a container of a height H along the z direction 
and the wall area L''"^, the derivative becomes the difference equation, and thus the discrete 
expression of the pressure is given by: 



where Z[H) is the partition function with the height H and Z{H — 1) is that with H-1. 
Note that Z{H — 1) is equivalent to Z{H) with an inflnite repulsive potential at the wall, in 
which case the particles cannot appoarch the wall at z=H and thus with U the interaction 
potential between the wall and particles, the original partition function, Z{H, U — 00), will 
reduce to Z(H — 1, [/ = 0). In order to utilize this observation in designing the Monte Carlo 
method, we introduce a parameter A = exp{—U/kT) with k the Boltzmann constant. Note 



and = N/V, we can easily obtain eq.(2). Note that the pressure of the lattice gas exhibits 




k{T/L'^-^)ln{Z{H)/Z{H - 1)) 



(6) 
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that A = corresponds to U — oo, and X — 1 to U — 0. Let the total number of particles 
in the system be N, and Nyj be the particles at the wall. The partition function, Z, for this 
modified system becomes: 

Z{H, X) = J2 exp{-E/kT) = A^»0(Ar^, N - N^) (7) 

E N^=0 

where the summation is taken over all possible energy configurations and ft is the total 
number of configurations with A^^ particles at the wall and N — N^^ particles in the bulk. 
The probablity, P{N^), of finding particles at the wall in the presence of the potential 
U is, P{N^) = X^-'fl{N^, N - N^)/Z. Further, it is obvious that Z{H) = Z{H, A = 1) and 
Z{H — 1) — Z{H,X — 0). Therefore, if the mean number of particles at the wall is denoted 
by < Nw >, then it is given by: 

, dlnZ 

< >= (8) 
Next, we use a simple mathematical identity to express the ratio ln{Z[H)/Z[H — 1)) in 
terms of < Ny, > as: 

ln(Z(H) iZiH - 1)) = IniZiH, A = 1) iZiH, A = 0)) = /' rfA^^^^^J^^ (9) 

Jo aX 

Hence, we obtain the desired expression for the pressure^: 

P/kT^['^p^ (10) 
Jo A 

where =< > / L'^~^ is the density at the wall. 

We now demonstrate here that Eq.(lO) is indeed identical to Eqs.(3) and (5). Consider 
first the case where the excluded volume interactions are present. Since the partition function 
Z (Eq.(7)) contains only the configurational term, it is straightforward to write down the 
expression for Z for particles in a volume V = HA = H • L^'^~^^ with A = L'^'"^ the area 
of the container. We first consider the case with excluded volume. The partition function 
becomes, 

Z{H,X)= X''-n{A,N^)n{B,N -N^) (11) 

Ny.=0 



where B = {H — 1)L('^~^). The first term,A^"' in (A-1) is due to the wall potential, and the 
second and the third terms are the total number of ways of putting Ny^ particles in the wall 
of area A and N — Nyj particles in the bulk of a volume B. Both Q's in the summation 
are given by the binomial coefficient, i.e.namely Q{Q, R) — Q\/R\{Q — R)\. Let the density 
of the particle, N/B = < 1. Then, in the thermodynamic limit, Nyj « N. Using the 
Sterling's formula, we find: 



lnn{B, N-N„)^ -B[(t>ln(t> + (1 - 0)Zn(l - 0)] + Njn[(t>/{1 - 0)] + 0{N^/B) 
Hence, 

N-N^)^ exp{Bso) [0/(1 - (12) 

where So is the bulk entropy per site; i.e. So — —[(plrKf) + (1 — (f))ln{l — 0)]. By putting 
(12) into (11), we now obtain the closed expression for the partition function Z: 

Z = exp(Ss„)[l + A0/(l-0)]^ 
from which we obtain the exact formula for the density at the wall, pw'. 

= {X/A)^ = Aa(0)/(1 + Xam (13) 

with a{4>) = 0/(1 — 0). Note that Pu,(A = 1) = as expected. Hence, we find 

P/kT = dXp^/X = -ln(l - 0) (14) 

In the case when the excluded volume interaction is not present, the partition function 
becomes much simpler: 

N 

Z= J2 X^'-niN, Ar^)A^-S^-^- = S^(l + XA/B)"" (15) 
Hence, in the thermodynamic limit of ^4/5 << 1, the density at the wall, pw, becomes: 
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= {\/A) 



dlnZ 
dX 



1 + XA/B 



X(j) 



Xcf) 



(16) 



from which follows the ideal gas law: 



PjkT= [ 



1 dX 



(17) 



J o 



IV. Random Walk Simulations 

We now carry out Monte Cairo random walk simulations to measure the density at the 
wall and compare them with the formulas (13) and (16). 

We take a two dimensional lattice of size H — 20 and L — 10 with the volume V — HL. 
We now turn on the repulsive wall potential < U < oo at z—R, and introduce a parameter 
< A = exp{—U/T) < 1. The wall potential is short range and thus acts only when 
the particle is right next to the wall, i.e., a,t z — H — 1. At z — H — 1, the particle 
moves to the wall ai z — H with the probability p — X. This is the standard metropolitan 
algorithm. Otherwise, the particles are free to move in the bulk (with the exception of 
the excluded volume interaction). We also impose periodic boundary conditions along the 
vertical axis. With these simple rules, the Monte Carlo simulations proceed as follows. 
Initially N particles with the bulk density 4> = N/V are randomly placed on a lattice and 
they undergo random walk. With the excluded volume interaction, multiple occupancy on 
a lattice point is prohibited. Without the excluded volume interaction, multiple occupancy 
is allowed. When a particle is right next to the right wall, then it moves there with the 
probability A. When an attempt has been made to move all the N particles once, then one 
Monte Carlo(MC) step is elapsed. At each MC step, the number of particles at z — H — 1 
is registered. Then, after M MC time steps, the density of the particles near the wall, 
Pu; Bi z — H — 1 is computed as the average over MC time t and configurations, namely: 
Pw — [J2t Ny,{t)/M]/L with Ny,{t) the muber of particles at z — H — 1 at b, given MC time 
t. is then measured as a function of A for different bulk density (p's. We now present 
simulation data. 
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We measured the particle density at the wall, pyj{\,(f)), with and without excluded in- 
teractions, for the bulk density = 0.1, 0.3, 0.5, 0.7, 0.9 and for A = 0.1, 0.3, 0.5, 0.7, 0.9, 1.0. 
By symmetry, pwaiii^i 4>) = 4>- The simulation data are presented in Fig.l for the excluded 
volume interaction and in Fig.2 without the excluded volume interactions. The dotted lines 
in both Figures are the predicted formulas, (13) and (16), and the squares are the simulation 
data. Even for a small size of lattice used in the simulations, the agreement between simu- 
lation data and the predicted forumla appear to be quite remarkable. The typical error bars 
are less than a few percent and the error bars are expected to vanish in the thermodynamic 
limit as the size increases. The good agreement between the simulations and the prediction, 
even for a small size of lattice, is precisely because of the fact that the ratio A/B, in equa- 
tion (16) is essentially the ratio of one dimensional lattice sites over those of two dimensions, 
which is in the presence case of order 1/20 = 0.05 << 1. Even for a modest size of two 
dimensional lattice, this ratio is small and its contribution is almost insignificant in higher 
dimensions. Further, the logarithmic singularity in the pressure at the closed packed density 
is also well captured in the lattice simulations with the typical error bars are less than a few 
percent. We now conclude with two comments. 

First, while this method appears to produce excellent results, the method becomes a 
little problematic when there is an external bias, say the gravity toward the wall. In this 
case, unless the bias is small, it will eventually pull all the particles to the wall, giving the 
wall density essentially close to one, and largely independent of A. In this case, it becomes 
quite difficult to exploit the relation (10) for diffent A's. We have not come up with a method 
to overcome this difficulty. Such a method, if successful, will be quite useful because it will 
enable one to directly measure the force profile of granular materials in a container, which 
is the subject of current studies [9-12]. 

Second, the method, however, is expected to work quite well in the absence of bias for 
other interacting systems. An example would be the charged or uncharged lattice gas with 
either long or short range interactions such as Coulomb interactions, or other Ising type 
interactions near the critical point. Such studies will require elaborate finite size scaling 
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analysis and extensive large scale simulations, which will be reported in future communica- 
tions. 
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Figure captions: 



Fig. 1. Comparison between the Monte Carlo data and the exact formula (13) for the lattice 
gas with the excluded volume interaction for given bulk density 0's for different X's. The 
vertical axis is the particle density at the wall, pw, and the horizontal axis is A. From 
the bottom to top, (f) = 0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9. Dotted lines are the exact 
formula(13). 

Fig. 2. Same as in Fig.l except that multiple occupancy is allowed at the lattice sites. 
Dotted lines are the exact formula(16). 
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